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1  ARO  Report 

1.1  Modelling  of  Parachute  Canopy 

Simulation  of  parachute  inflation  via  computational  methods  has  attracted  the  attention  of 
scientists  at  the  Department  of  Defense  laboratories  and  academia  alike.  Some  of  the  suc¬ 
cessful  studies  include  those  reported  in  [28,  27,  24,  29,  25,  26],  where  the  Deforming-Spatial- 
Domain/Stabilized  Space-Time  (DSD/SST)  method  [34,  35]  was  used  in  [27,  24,  29,  25,  26]. 
The  studies  in  [37,  36,  31,  32,  38]  also  used  the  DSD/SST  as  the  core  numerical  method,  but 
involved  new  versions  and  special  techniques.  These  studies  successfully  addressed  the  com¬ 
putational  challenges  in  handling  the  geometric  complexities  of  the  parachute  canopy  and  the 
contact  between  parachutes  in  a  cluster.  Kim  and  Peskin  et  al.  used  the  immersed  bound¬ 
ary  method  to  study  the  semi-opened  parachute  in  both  two  and  three  dimensions  [10,  12], 
their  simulations  are  on  small  Reynolds  number  (about  300)  and  applied  to  parachute  with 
payload  of  several  grams.  Karagiozis  used  the  large-eddie  simulation  to  study  the  parachute 
in  Mach  2  supersonic  flow  [7].  Purvis  [20,  21]  used  springs  to  represent  the  structures  of  the 
fore-body,  the  suspension  lines,  and  the  canopy,  etc.  In  these  papers,  the  authors  used  cylin¬ 
drical  coordinates  with  the  center  line  as  the  axis.  In  the  paper  by  Strickland  et  al.  [30],  the 
authors  developed  an  algorithm  called  PURL  to  couple  the  structure  dynamics  (PRESTO) 
and  fluid  mechanics  (CURL),  in  which  mass  is  added  to  each  of  the  structure  node  based 
on  the  diagonally  added  mass  matrix  and  a  pseudo  is  computed  from  the  fluid  code  which  is 
the  sum  of  the  actual  pressure  and  the  pressure  associated  with  the  diagonally  added  mass. 
Tutt  and  Taylor  [41,  40]  simulated  the  parachute  through  the  LS-DYNA  code.  They  used  an 
Eulerian-Lagrangian  penalty  coupling  algorithm  and  multi-material  ALE  capabilities  with 
LS-DYNA  to  replicate  the  inflation  of  small  round  canopies  in  a  water  tunnel. 

Our  computational  approach  to  the  simulation  of  the  parachute  delivery  system  is  based 
on  the  front  tracking  platform.  In  the  last  ARO  funded  project  during  which  the  PI  served 
as  the  principal  developer,  we  established  a  computational  platform  for  the  parachute  study 
by  employing  the  spring  model  for  the  parachute  canopy  and  the  string  chords.  We  designed 
a  set  of  new  data  structure  to  allow  the  application  of  the  FronTier  library  to  track  the 
dynamic  motion  of  fabric  surface  driven  by  the  gravitational  force  of  payload  and  the  fluid 
pressure.  We  discretize  the  fabric  surface  into  a  homogeneously  triangulated  surface  mesh. 

1.1.1  Analysis  of  the  Spring  System 

We  made  a  detailed  analysis  of  the  spring  system  which  is  used  to  model  the  fabric  surface 
of  the  parachute  canopy  (Figure  1).  When  no  external  driving  force  is  applied,  the  fabric 
surface,  which  is  represented  by  the  spring-mass  system,  is  a  conservative  system  whose 
total  energy  (kinetic  energy  plus  potential  energy)  is  a  constant.  Assuming  each  mesh  point 
represents  a  point  mass  m  in  the  spring  system  with  position  vector  X;,  the  kinetic  energy 
of  the  mass  point  i  is  Tj  =  |m|xj|2,  where  x,  is  the  time  derivative,  or  velocity  vector  of  the 
mass  point  i.  Here  we  report  the  work  under  the  support  of  ARO  grant  W911NF1310291  on 
the  improvement  of  the  spring  model,  its  verification  of  Young’s  modulus  and  Poisson  ratio, 
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Figure  1:  The  spring  model  on  a  triangulated  mesh.  Each  vertex  point  in  the  mesh  represents 
a  mass  point  with  point  mass  m.  Each  edge  of  the  triangles  has  an  equilibrium  length  set 
during  initialization  and  the  changing  length  exerts  a  spring  force  on  the  two  neighboring 
vertices  in  opposite  directions.  Gores  are  added  as  curves  with  larger  spring  constant.  This 
plot  shows  the  mesh  of  canopy  with  gores. 

and  the  validation  in  the  simulation  of  parachute  inflation  and  descending. 

The  original  spring  model  [17]  that  serves  as  the  basis  for  this  effort  is  a  simplified 
mesoscale  model  which  assumes  the  force  required  to  bend  the  surface  is  negligible  and  the 
force  to  stretch  the  surface  is  proportional  to  the  displacement  from  the  equilibrium  distance 
between  adjacent  mass  points.  In  this  model,  the  kinetic  energy  and  the  potential  energy  of 
the  triangulated  mesh  is  given  by 

1  1  N  N 

Ti=-m,|X,|2,  V=j^EidXi  (1) 

i=  1  j= 1 

where  X,  is  the  position  of  the  vertex  i,  k  is  the  spring  constant,  is  the  equilibrium 
length  of  the  side  shared  by  vertices  X*  and  Xj,  and  rjVJ  is  a  boolean  variable  for  adjacency. 
Delingette  modification  is  added  to  the  spring  model.  This  has  made  possible  for  a  direct 
link  between  the  meso-scale  spring  model  and  the  elastic  membrane  model  in  continuum 
mechanics. 

Illustrative  features  of  Delingette’s  model  are  given  in  Figure  2.  The  energy  of  membrane 
W(Tx q)  that  is  required  to  deform  a  single  triangle  TXo  with  vertices  {X10,  X2o,  X30}  into 
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its  deformed  position  Tx  with  vertices  {X^,  X2,  X3}  consists  of  two  parts, 

•  The  energy  of  three  tensile  springs  that  prevent  edges  from  stretching. 

•  The  energy  of  three  angular  springs  that  prevent  any  change  of  vertex  angles. 

For  a  triangle  in  equilibrium  Tx0,  the  initial  states  are  given  by  area  ^4Xo ,  angles  cq,  and 
lengths  If  (j  6  1,2,  3)  in  equilibrium  while  Ax,  A  and  denote  the  area,  angles,  and  lengths 
of  the  deformed  triangle  Tx,  respectively. 

The  edge  elongation  can  be  written  as  dli  —  lt  —  l °.  The  potential  energy  is  given  [2]  as 

3  1  3 

w(t*„)  =  y.  Yx°  w2  +  X 

i= 1  Z  i= 1 

j=(i+ 1)  (mod  3) 
k=(i+ 2)  (mod  3) 


where 


is  the  tensile  stiffness  and 


(7°)2(2  cot2  cq( A  +  fi)  +  fi) 

8Ax„ 


rx  Ijlki 2  c°t  c°t  +  /^)  —  aO 

8Ax0 


(2) 


is  the  angular  stiffness,  where  j  —  (i  +  1)  (mod  3)  and  k  —  (i  +  2)  (mod  3).  7  and  /i  are  the 
Lame  coefficients  of  the  material.  These  coefficients  are  simply  related  to  the  two  physically 
meaningful  parameters  defined  in  planar  elasticity  for  a  membrane,  that  is,  Young’s  modulus 
E  and  the  Poisson  ratio  v  [3]: 


A 


Ev 
1  —  v2 


and  /i 


E{  1  -  is) 
1  —  is2 


Young’s  modulus  quantihes  the  stiffness  of  the  material,  whereas  the  Poisson  ratio  charac¬ 
terizes  the  material  compressibility. 

Through  the  application  of  Raylcigh-Ritz  analysis  the  fabric  surface,  represented  by  the 
triangular  mesh,  should  evolve  by  minimizing  its  membrane  energy.  Therefore,  along  the 
opposite  derivative  of  that  energy  with  respect  to  the  nodes  of  the  system,  that  is,  the 
deformed  positions  Xp 


The  membrane  deformation  energy  of  the  whole  triangulation  is  the  sum  of  the  energy  of 
each  triangle.  Thus,  we  obtain  the  force  at  each  vertex  point  as  Eq.  (4). 
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Figure  2:  (a)  Rest  triangle  7x0  whose  vertices  are  Xi0.  (b)  Deformed  triangle  7x  whose 
vertices  are  X,;. 


Figure  3:  Triangles  Xf  and  T2  share  X,  and  Xj,  the  other  vertices  of  triangles  T\  and  T2  are 
Xm  and  Xn  respectively. 


N 


Fj  'y^j  VijFjj 

3= 1 


(4) 


As  illustrated  by  Figure  3,  if  X,;  and  Xj  are  shared  by  two  triangles  Tf  and  T2,  and  the 
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other  vertices  of  triangles  Ti  and  T2  as  denoted  as  Xm  and  Xn  respectively,  we  have 
Fjj  =  (( kjj  +  *£  )dlij  +  (' y^dlim  +  'y^dljm  +  7 J2dhn  +  7 'J2(Mjn))eij 

kijdlijGij  'YijdlijGij  (5) 

where  k^  =  kjf  +  kj?,  =  (yf1  dlim  +  rfldljm  +  7 J2dlin  +  7 J2dljn)/dlij  and  is  the  unit 
vector  from  X?:  to  Xj. 

If  the  second  term  in  Eq.  (5)  is  neglected,  Dclingette’s  model  is  the  same  as  model  used 
in  [9]  except  that  the  spring  constant  varies  if  the  corresponding  initial  triangles  deviate 
from  isosceles  triangles.  Numerical  evidence  suggests  that  both  the  variation  of  k^  and  the 
modification  from  the  second  term  (due  to  angular  stiffness)  are  significant  in  the  above 
cases. 


Figure  4:  Convergence  test  of  the  drum  membrane  under  mesh  refinement.  From  left  to  right 
the  computational  mesh  of  the  domain  are  15,  30,  and  60  respectively.  The  total  mass  of  the 
membrane  is  kept  constant  in  the  simulations.  The  upper  three  plots  show  the  membrane 
position  at  t  —  1  and  the  lower  plots  show  the  membrane  position  at  t  —  2. 


1.1.2  Numerical  Convergence  Verification 

Numerical  convergence  under  computational  mesh  refinement  is  a  crucial  step  in  assessing 
the  validity  of  spring-mass  model.  For  this  purpose,  we  carried  out  two  sequences  of  simu- 
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reference  mesh  size 

eA 

Cp 

15  and  30 

0.02528 

1.91810 

1.97967 

30  and  60 

0.01507 

1.21784 

1.21772 

60  and  120 

0.00604 

0.53550 

0.52837 

Table  1:  Convergence  tests  of  spring  model  for  a  fabric  drum.  In  the  computational  se¬ 
quences,  the  total  mass  of  the  membrane  is  fixed.  As  the  number  of  points  increases,  the 
point  mass  is  reduced  accordingly.  Cauchy  error  is  calculated  on  two  consecutive  mesh  se¬ 
quences.  Column  6a,  e*,,  and  ep  are  errors  of  total  area,  total  kinetic  energy  and  total  spring 
potential  energy  respectively.  The  numerical  results  show  the  first  order  convergence. 


lations  with  increasing  number  of  grid  points.  Our  convergence  study  includes  both  a  one 
dimensional  string  in  the  two  dimensional  plane  and  a  two  dimensional  membrane  in  the 
three  dimensional  space. 

In  one  case  of  the  verification  study,  a  circular  vibrating  membrane  with  radius  of  0.4m 
and  total  weight  of  380 g  is  considered.  The  membrane  is  linearly  perturbed  at  initial  time 
and  its  circular  boundary  is  fixed.  A  sequence  of  four  cases  with  computational  mesh  153, 
303,  603,  and  1203  are  computed.  The  errors  of  total  area  are  shown  in  Table  1.  From 
Table  1,  we  can  clearly  see  that  the  Cauchy  errors  are  decreasing  as  the  computational  mesh 
is  refined.  The  convergence  rate  for  the  membrane  is  also  of  first  order. 


1.1.3  Verification  of  Young’s  Modulus  and  Poisson’s  Ratio 

Young’s  modulus,  also  known  as  the  tensile  modulus  or  elastic  modulus,  is  a  measure  of  the 
stiffness  of  an  elastic  material  and  is  a  parameter  used  to  characterize  materials.  It  is  defined 
as  the  ratio  of  the  stress  along  an  axis  over  the  strain  along  that  axis  in  the  range  of  stress 
in  which  Hooke’s  law  is  valid.  In  solid  mechanics,  the  slope  of  the  stress-strain  curve  at 
any  point  is  called  the  tangent  modulus.  The  tangent  modulus  of  the  initial,  linear  portion 
of  a  stress-strain  curve  is  called  Young’s  modulus.  Young’s  modulus,  E ,  can  be  calculated 
by  dividing  the  tensile  stress  by  the  tensile  strain  in  the  elastic  portion  of  the  stress-strain 
curve: 

E  tensile  stress  a  F/A0  Fl0 

tensile  strain  e  A l/l0  A0Al 

where  E  is  Young’s  modulus,  F  is  the  force  exerted  on  the  object  under  tension,  A0  is  the 
original  cross-sectional  area  through  which  the  force  is  applied,  A l  is  the  change  of  length  of 
the  object  from  the  original  length  /0  of  the  object. 

Poisson’s  Ratio  is  the  negative  ratio  of  transverse  strain  to  axial  strain.  When  a  material 
is  stretched,  it  usually  tends  to  contract  in  the  directions  transverse  to  the  direction  of 
stretching.  The  Poisson  ratio  is  the  ratio  of  relative  contraction  to  relative  stretching. 


de 


v  = 


trans 


de, 


axial 


dey 

dex 


(7) 
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where  v  is  the  Poisson’s  Ratio,  etrans  is  transverse  strain  and  eaxiai  is  axial  strain.  Theoret¬ 
ically,  in  the  case  of  small  deformations  Poisson’s  Ratio  v  was  computed  by  Eq.  (8)  and  in 
the  case  of  large  deformation  it  was  computed  by  Eq.  (9). 


Aw/ wq 

Al/lo 


(8) 


V  =  -Zo0(i+Aj/to)(l  +  a  w/w)  (9) 

To  verify  that  the  spring  model  can  catch  isotropic  elastic  material’s  Young  modulus  and 
Poisson  ratio,  we  carried  out  a  set  of  simulations  by  stretching  fabric  surfaces  with  different 
Young’s  modulus  and  Poisson  ratios.  These  simulations  start  with  a  fabric  surface  which 
has  its  original  length  l0  =  0.1  m  and  original  width  wo  =  0.02 m.  The  fabric  is  then  pulled 
along  the  direction  of  the  longer  side  of  it  with  a  distributed  force. 

Firstly,  three  groups  of  simulations  which  stretch  fabric  surfaces  with  Poisson  ratios  of 
—0.14,  —0.22,  and  —0.30,  respectively  are  carried  out.  Each  group  with  Eve  different  values 
of  Young’s  modulus,  that  is  0.1  GPa,  0.2 GPa,  0.3 GPa,  OAGPa,  and  0.5 GPa  respectively. 
The  fabric  surface  is  pulled  at  one  end.  Young’s  modulus  and  Poisson  ratio  are  calculated 
from  the  deformation  of  the  surface  and  the  force  added  on  it.  The  results  are  summarized 
by  Table  2  which  shows  that  the  spring-mass  model  nicely  reproduces  the  values  of  Young’s 
modulus  and  the  Poisson  ratio  from  the  input. 

Next,  a  group  of  simulations,  stretching  the  fabric  surface  whose  Young’s  modulus  and 
Poisson  ratio  are  fixed  at  0.5 Gpa  and  —0.14  respectively,  are  carried  out.  The  total  number 
of  triangles  of  these  triangulations  use  in  the  simulations  changes  from  1143,  1590,  2468, 
to  4520.  The  strain  of  the  elongation  change  from  0.002,  0.004,  0.006,  0.008  to  0.01.  The 
measured  Young’s  modulus  and  Poisson  ratio  from  these  simulations  are  demonstrated  in 
Figure  5. 

The  numerical  solutions  imply  that  the  revised  spring-mass  model  based  on  the  derivation 
by  Delingette  [2]  can  accurately  simulate  an  isotropic  elastic  membrane  in  the  linear  region 
with  strain  up  to  0.01.  The  error  between  the  spring  model  and  continuum  model  increases 
when  elongation  is  too  large  and  the  deformation  reaches  the  nonlinear  regime.  For  parachute 
simulation,  the  strain  of  the  canopy,  even  during  the  most  dynamic  phase  of  inflation,  should 
still  be  in  the  linear  regime.  Therefore  the  spring  model  is  an  excellent  model  for  such 
simulations. 

It  should  be  mentioned  that  in  the  case  in  which  the  fabric  surface  is  compressed  and 
it  adjusts  itself  with  wrinkles,  the  convergence  is  not  obvious.  In  some  cases,  a  small  per¬ 
turbation  of  the  initial  condition  may  result  in  substantially  different  folding  and  wrinkling 
patterns.  How  to  describe  such  case  and  define  its  mathematical  convergence  remain  as  an 
open  question  to  the  model. 


1.1.4  Fluid  Canopy  Coupling 

To  model  the  parachute  system,  an  accurate  coupling  between  the  Naiver-Stokes  equation 
and  the  structure  dynamics  must  be  carefully  considered  near  the  canopy  surface.  The 
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Group  1 

Young’s 

Modulus  results  ( GPa ) 

Poisson  Ratio  results 

Input 

Numerical 

Input 

Numerical 

case  1 

0.1 

0.0963016606 

-0.14 

-0.141887 

case  2 

0.2 

0.1926061554 

-0.14 

-0.141920 

case  3 

0.3 

0.2888656609 

-0.14 

-0.141910 

case  4 

0.4 

0.3852341348 

-0.14 

-0.141913 

case  5 

0.5 

0.4815205552 

-0.14 

-0.141905 

Group  2 

Young’s 

Modulus  results  (GPa) 

Poisson  Ratio  results 

Input 

Numerical 

Input 

Numerical 

case  1 

0.1 

0.0957361901 

-0.22 

-0.223342 

case  2 

0.2 

0.1914918931 

-0.22 

-0.223332 

case  3 

0.3 

0.2870650506 

-0.22 

-0.223382 

case  4 

0.4 

0.3829895385 

-0.22 

-0.223349 

case  5 

0.5 

0.4784196259 

-0.22 

-0.223365 

Group  3 

Young’s 

Modulus  results  (GPa) 

Poisson  Ratio  results 

Input 

Numerical 

Input 

Numerical 

case  1 

0.1 

0.0953452395 

-0.30 

-0.305362 

case  2 

0.2 

0.1906619377 

-0.30 

-0.305342 

case  3 

0.3 

0.2860538596 

-0.30 

-0.305343 

case  4 

0.4 

0.3814011184 

-0.30 

-0.305359 

case  5 

0.5 

0.4767297379 

-0.30 

-0.305355 

Tabic  2:  Young’s  modulus  and  Poisson  ratio  verification 


Figure  5:  Young’s  Modulus  (Left)  and  Poisson’s  Ratio  (Right).  We  tested  the  spring  model 
by  stretching  the  fabric  surface  to  different  lengths.  The  numerical  results  show  that  the 
spring-mass  model  catches  the  fabric’s  Young’s  Modulus  and  Poisson  ratio  nicely  in  the  linear 
regime  of  strain. 


method  we  designed  for  the  simulations  of  air  delivery  system  uses  the  superposition  of 
impulse  on  every  mass  point.  Each  mass  point  in  the  spring  system  acts  as  an  clastic 
boundary  point  and  exerts  an  impulse  on  the  fluid  in  its  normal  direction.  Our  algorithm 
ensures  that  the  action  and  reaction  between  the  spring-mass  point  and  the  fluid  solver  are 
equal  in  magnitude  and  opposite  in  directions,  a  requirement  of  Newton’s  third  law. 

We  applied  Peskin’s  delta  function  immersed  boundary  method  [10,  12] 

f(x,  t)  —  I  F(s,  t)5(x.  —  X(s,  t))ds.  (10) 

The  difference  between  our  method  and  Kim  and  Peskin’s  method  lies  in  the  calculation  of 
F.  Instead  of  computing  the  tension  through  the  derivative  with  respect  to  the  arc  length, 
we  use  the  impulse  of  the  mass  point  as  a  result  of  the  superposition  of  three  forces  from  the 
spring  system,  that  is 

F(xi,i)  =  d(Ig  +  Ip  +  Is)/dt  (11) 

Eq.  (11)  is  more  physically  realistic,  especially  because  ls  is  obtained  from  the  spring  equa¬ 
tions.  In  the  canopy  spring  system,  we  have  observed  that  the  tension  is  high  at  the  top  of 
the  canopy  where  the  curvature  is  almost  zero. 

1.2  The  k  —  £  Turbulence  Model 

A  round  parachute  canopy  poses  several  challenging  problems  for  the  understanding  of  bluff- 
body  wake,  which  has  been  studied  extensively  during  the  past  decades  [18,  19,  33,  8]. 
However,  the  main  difference  between  fabric  canopies  and  rigid  bluff  bodies  are  fabric  flex¬ 
ibility  and  permeability  [5].  Thus  the  coupling  of  the  vertical  flow  with  the  permeable 
and  flexible  fabric  may  have  significant  difference  with  rigid  bluff  bodies.  Although  direct 
numerical  simulation  (DNS)  and  large  eddy  simulation  (LES)  may  give  a  relatively  accu¬ 
rate  prediction  of  turbulence  behavior,  the  high  computational  cost  are  still  obstacles  of 
extending  these  methods  to  high  Reynolds  number,  which  typically  exceeds  several  millions 
during  parachute  inflation.  Therefore,  eddy  viscosity  models  based  on  the  Reynolds  Aver¬ 
aged  Navier-Stokes  (RANS)  equations  are  commonly  used  in  turbulence  simulation.  One  of 
the  most  popular  RANS-based  turbulence  model  is  the  k  —  £  family,  which  automatically 
provide  the  turbulence  length  scale  or  its  equivalent  and  are  thus  complete  [43].  This  kind 
of  model  still  involved  some  empirical  constants,  which  are  widely  used  and  fitted  well  with 
some  experimental  data. 

In  this  section,  we  firstly  give  a  brief  introduction  to  the  k  —  £  family,  including  stan¬ 
dard,  RNG  and  realizable  model.  Then  a  series  of  simulations  were  performed  to  make  a 
comparison  between  these  models. 

1.2.1  Mathematical  Formulation 

Following  the  Boussinesq  assumption  [43] ,  the  hydrodynamic  behavior  of  a  turbulent  incom¬ 
pressible  fluid  is  governed  by  the  RANS  equations  for  the  mean  velocity  U  and  pressure 
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p 


dU 

~dt 


+  u  •  w 


-Vp  +  V  •  {{u  +  uT)[VU  +  Vf/T]) 

V-U  =  0 


(12) 


where  //  is  the  kinematic  viscosity  and  vt  is  the  turbulent  eddy  viscosity  which  is  supposed 
to  approximate  the  turbulent  fluctuations.  In  standard  k  —  e  model,  the  eddy  viscosity  is 
defined  as 

vt  =  — ,  (13) 

where  k  is  the  turbulence  kinetic  energy  and  e  is  the  dissipation  rate.  To  compute  k  and  e, 
two  additional  convection-diffusion-reaction  equations  are  defined  as: 

^  +  V-{kU-{u+  V-f)Wk)  =Pk-e  (14) 

at  ok 

^  +  V  -{eU-  (u  +  Uf)S7e)  =  |(C'1Pfc  -  C2e)  (15) 

where  Pk  =  ^|V£/  +  V UT | 2  is  the  production  of  turbulent  kinetic  energy.  For  the  standard 
k  —  e  model,  the  default  values  of  the  involved  empirical  constants  are:  C ^  =  0.09,  C\  = 
1.44,  C2  =  1.92,  8k  =  1.0,  5e  =  1.3.  Although  simple  and  efficient,  the  standard  model  is 
unable  to  capture  the  effects  of  smaller  scales  of  motion  due  to  its  single  turbulence  length 
scale.  In  order  to  account  for  the  different  scales  of  motion,  a  mathematical  technique,  Re- 
Normalization  Group  (RNG)  method,  was  used  to  derive  a  turbulence  model  similar  to  the 
standard  one,  results  in  a  modified  form  of  the  £  equation: 


^  +  V  •(£[/-  (u  +  Vf)S/£)  =  |(C1Pfc  -  C*2£)  (16) 


C*  =  C2  + 


~  v/vo) 

1  +  f3r)3 


(17) 


rj  =  kS/£,  S  is  modulus  of  the  mean  rate  of  strain  tensor.  The  coefficients  are  derived 
explicitly  in  the  RNG  procedure  and  are  also  listed  here  for  completeness:  C =  0.0845, 
Ci  =  1.42,  C2  =  1.68,  Sk  =  0.7194,  5e  =  0.7194  [44], 

Using  standard  or  RNG  model  may  produce  unphysical  results  when  the  mean  strain 
rate  is  large.  The  most  straightforward  way  to  ensure  the  realizability  is  to  make  C ^  variable 
according  to  the  mean  flow.  It  is  around  0.05  in  homogeneous  shear  flow  and  around  0.09  in 
the  inertial  sub-layer  of  boundary  layer  flow.  The  kinetic  energy  equation  remains  the  same 
as  standard  model,  while  the  equation  for  dissipation  rate  is  modified  as  follows: 


%  +  V  •  (eU  -  (v  +  ^)Ve)  =  UC.S  -  C2 


dt 


5f 


k 


k  +  y/[ve) ' 


(18) 


77 

max{0.43, - r? 

77+5 


£ 


(19) 
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The  coefficient  C),  is  determined  to  meet  the  positive  normal  stress  and  Schwarz’s  inequalities. 
The  details  of  implementing  Realizable  model  can  be  found  in  [23]  and  the  reference  therein. 

The  implementation  of  k  —  £  model  has  no  difference  with  solving  other  convection- 
diffusion-reaction  equations,  except  that  two  things  need  to  be  noticed.  One  must  specify  an 
appropriate  boundary  condition  to  a  solid  wall  for  the  velocity  and  two  turbulence  parame¬ 
ters.  Due  to  the  extremely  high  Reynolds  number  in  parachute  simulation,  it  is  worthwhile 
to  use  the  wall  functions  to  bridge  the  viscosity-affected  region  and  the  fully-turbulent  region 
avoiding  the  need  for  resolution  of  strong  velocity  gradients  [15].  Another  task  is  to  avoid 
loss  of  positivity  of  k  and  £  due  to  computational  errors.  This  can  be  achieved  by  simply 
keeping  the  coefficients  positive  for  the  linearized  equations  without  touching  any  primitive 
variables  [16].  In  the  next  subsection,  a  series  of  parachute  simulation  are  performed  to 
compare  the  results  between  these  three  different  models. 

1.2.2  Numerical  Implementation  and  Test  Results 

In  this  section,  we  examine  incompressible  flow  past  a  canopy  surface  placed  in  a  channel 
orthogonal  to  the  streamwise  fluid.  The  numerical  domain  is  set  to  be  12m  x  24m  with 
Dirichlet  boundary  conditions  in  the  streamwise  while  slip,  no  penetration  boundary  condi¬ 
tions  are  applied  at  the  channel  walls  and  canopy  surfaces.  The  configurations  are  illustrated 
in  Figure  6.  The  Reynolds  number  is  defined  as: 

V 

where  U  =  3ms -1  is  the  inflow  velocity,  L  =  7.0m  is  the  width  of  the  canopy,  v  —  1.5  x 
10^5m2s^1  is  the  kinematic  viscosity  of  the  air,  thus  Re  ~  106. 

We  use  finite  difference  method  to  solve  the  RANS  equations.  The  pressure  and  velocity 
in  equation  (12)  are  decoupled  with  Chorin’s  projection  method  and  advanced  in  small 
time  step  according  to  4-th  order  Runge-Kutta  method  with  numerical  flux  evaluated  by 
5-th  order  WENO  scheme.  The  equations  for  k  and  £  are  discretized  with  Crank-Nicolson 
scheme  to  achieve  second  order  of  accuracy.  Some  numerical  techniques  similar  to  [16]  are 
applied  to  preserve  the  positivity  of  turbulence  parameters. 

We  show  the  numerical  results  by  applying  different  turbulence  models.  Figure  7  and 
Figure  8  display  the  viscosity  and  velocity  computed  from  different  models.  The  RNG  and 
realizable  model  perform  better  than  the  standard  one  since  they  successfully  capture  the 
separation  and  rotation  of  the  air  flow. 

1.3  Numerical  Implementation  of  Porosity 

Fabric  porosity  of  parachute  canopy  greatly  affects  parachute  performance.  It  has  long  been 
known  that  fabric  permeability  is  an  important  factor  in  parachute  design.  By  carefully  con¬ 
sidering  the  balance  between  payload  and  drag,  a  finite  permeability  will  make  the  parachute 
much  more  stable.  The  permeability  of  the  parachute  material  often  plays  a  vital  role  in 
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12m 


Figure  6:  Experimental  settings 

the  design  of  parachute.  A  state-of-art  tuning  between  an  impervious  material  to  a  highly 
permeable  fabric  can  make  a  parachute  from  a  wandering  sloth  into  a  plummeting  stabilizer. 

There  are  two  forms  of  parachute  porosity:  for  canopy  manufactured  from  solid  fabric, 
the  porosity  (sometimes  referred  as  fabric  permeability)  is  defined  as  the  airflow  through 
canopy  cloth  in  ft3/ ft2 /min  (cubic  feet  per  minute  per  square  feet)  at  /  inch  water  pres¬ 
sure;  for  slotted  canopies  which  has  geometric  openings,  porosity  is  defined  in  percent  as 
the  ratio  of  all  open  areas  to  the  total  canopy  area.  Most  personnel  parachutes  and  the 
main  descent  parachutes  for  air  vehicles  use  materials  with  porosity  values  from  80  to  150 
ft? / ft2 /min.  Gliding  parachutes  use  almost  imporous  materials  from  0  to  5  ft? / ft2 /min. 
Slotted  parachutes  use  geometric  porosity  values  in  10%  to  35%  range  [14], 
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Figure  7:  Eddy  viscosity  from  different  models  at  t  —  10s,  from  left  to  right  are  standard, 
RNG  and  realizable  k  —  e  model 

Porosity  influences  parachute  drag,  stability  and  opening  forces.  Higher  porosity  de¬ 
creases  opening  forces  and  oscillation,  but  reduces  drag  forces  at  the  same  time  which  is 
usually  not  desirable.  As  a  result,  parachute  canopy  porosity  is  an  important  factor  to  con¬ 
sider  for  both  parachute  design  and  parachute  simulation.  Different  types  of  parachutes  have 
different  design  characteristics  requirements  [4,  14,  22], 

Given  the  significance  of  porosity,  there  are  several  attempts  to  simulate  porous  canopy 
parachute  motion.  Kim  and  Peskin  [11]  use  immersed  boundary  method  to  simulate  parachute 
motion  and  derive  the  relative  velocity  between  fluid  and  canopy  interface  using  Darcy’s  Law. 
Tutt  [39]  simulated  parachute  performance  under  LS-DYNA  and  use  Ergun  equation  to  de¬ 
scribes  the  magnitude  of  porous  flow.  Wang,  Aquelet,  Tutt,  Do,  Chen  and  Souli  [42]  simulate 
the  interaction  between  the  fluid  and  porous  medium  by  a  Euler- Lagrange  coupling  under 
LS-DYNA  framework. 

For  Kim  and  Peskin’s  method,  porosity  is  considered  as  the  leaking  pores  in  the  immersed 
elastic  boundary.  Let  fi  be  the  density  of  pores,  which  means  there  are  fids  pores  in  the 
interval  (s,  s  +  ds )  along  the  surface.  If  each  pore  has  an  aerodynamic  conductance  equal  to 
7,  which  means  the  flux  through  the  pore  is  y(pi  —  P2)  where  p\  and  P2  are  the  pressures  on 
the  two  sides  of  the  boundary,  then  the  flux  through  the  interval  (s,  s  +  ds)  of  the  boundary 
is  given  by  fi^(pi  —  P2)ds. 

Tutt’s  algorithm  is  closest  to  the  front  tracking  implementation.  Tutt  used  the  Ergun 
equation  to  describe  the  magnitude  of  porous  flow  velocity  at  a  given  pressure  difference 


13 


Figure  8:  Velocity  field  from  different  models  at  t  —  10s,  from  left  to  right  are  standard, 
RNG  and  realizable  k  —  £  model  with  color  indicating  the  magnitude  of  velocity 


based  on  two  coefficients  in  the  following  equation: 


A  P 
~L~ 


u  p 

K\  f  ^  Ii2 


v2  —  a  ■  Vf  +  b  ■  v2 


where: 


Ah 

I<2 


£ 3  ■  D2 

150  ■  (1  -e)2 

£3  ■  D 

1.75-  (1  -e) 


are  referred  to  as  the  viscous  and  inertial  factors  respectively.  D  is  defined  as  the  character¬ 
istic  length,  e  is  the  porosity  and  is  equal  to  the  ratio  of  the  void  and  total  volume.  Vf,  fi 
and  p  are  fluid  velocity,  viscosity  and  density  respectively.  In  our  FronTier-airfoil  code,  the 
pressure  difference  is  computed  after  projecting  the  velocity  into  the  divergence-free  space. 
We  can  then  use  the  pressure  difference  on  two  sides  of  the  canopy  to  compute  the  Vf  and 
add  it  to  the  advection  solver  as  flux  from  the  immersed  elastic  surface. 

A  cloth  porosity  of  27.4  ft3/ft2/min  at  |  inch  water  pressure  is  equivalent  to  1%  ge¬ 
ometric  porosity  [14].  This  enables  comparison  between  different  parachutes  and  uniform 
simulations  among  parachutes  which  have  different  porosity  types. 

We  introduce  the  concept  of  penetration  ratio  7,  which  is  a  dimensionless  coefficient 
0  <  7  <  1.  When  7  =  0,  it  means  no  penetration  will  happen  and  the  fluid  on  two  sides 
of  boundary  has  no  connection;  when  7  =  1,  it  means  the  interface  does  not  exit;  when 
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0  <  7  <  I,  part  of  the  fluid  can  penetrate  the  canopy  and  canopy  interface  is  a  porous 
medium.  The  fluid  held  is  described  by  Navier-Stokes  equation  and  solved  by  projection 
method  [1].  Coupling  penetration  ratio  into  projection  method  transforms  fluid  through 
porous  medium  simulation  into  boundary  condition  treatment. 

The  projection  method  contains  two  major  steps:  transportation  process  to  calculate 
intermediate  velocity  and  projection  process  to  derive  pressure  and  new  velocity. 


In  the  advection  part,  fifth-order  WENO  scheme  is  used  to  proceed  calculation  which 
requires  information  from  seven  nearby  points,  i.e. 


“h'  =  /(«?- 3XXlX« 


i+liUi+2)  ui+ 3)  • 


When  came  across  boundary,  interpolation  is  required  to  obtain  ghost  points.  Assume 
parachute  canopy  is  between  iq  and  ui+ 1,  then  uf+1,uf+2,u^+3  should  be  interpolated  as 

u^1host  =  9l({ul.J), 

U$°St  =  92  (R.,n})  , 

=  93  ( R.,n})  • 

Thus,  the  calculation  of  velocity  on  the  boundary  is 

„.k+l  _  f  („.k  k  k  k  k, ghost  k, ghost  k,ghost\ 

ui  —  J  \ui- 35  ui- 2?  ui- 1)  ui  5  ai+ 1  )  ui+2  >  ui+ 3  J  • 


Coupling  penetration  ratio  7,  we  set  boundary  points  as 

= a  -  7X1 + 

•*;r=a-7)4+TC" 
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h+T°  =  (1  —  7)  ui+3  +  7ui+3 


which  leads  to  new  construction  formula 


k+l  _  r  (  k  k  k  k  k,poro  k,poro  k,poro\ 

ui  —  J  \ui- 3?  ui- 2)  ai- 1?  ui  )  ui+l  ?  uj+2  )  ui+3  I  ' 

Diffusion  part  is  solved  using  Crank-Nicolson  scheme  which  has  common  form 


<+l  -  ±  1  /  -  2<+1  +  ng  <1  ~  2<  +  <-1 

At  2  1  Ax2  Ax2 


and  given  a  =  it  can  be  simplified  as 


— cmf+i  +  (1  +  2a)  u, •  —  au  Ax  =  au£+1  +  (1  —  2a)  uf  +  au";_1. 

Assume  there  is  canopy  interface  between  tp  and  Uj+i,  and  consider  no  porosity  case,  we 
get  the  velocity  from  interface  ukift^te .  Coupling  with  constant  velocity  boundary  condition, 
we  have 


k, state  1  o  \  k-\-l  /c+1  k, state  ,  /-i  o  \  /c  i  /c 

— a«Ai  +  (1  +  2a)  -  au^ x  =  aui|1  +  (1  -  2a)  u{  +  atq^ 


fc+l  k, state 


Introducing  penetration  ratio  gives 


-iau*+l  +  (1  +  2a)  uf+1  -  auf+11  =  2  (1  -  7)  a^ate  +  7auf+1  +  (1  -  2a)  uf  +  aujlj. 


k, state 


After  obtaining  intermediate  velocity,  projection  step  is  carried  out  to  obtain  pressure 
and  new  velocities.  Solving  Poisson  equation  gives  pressure 

Pi+ 1  -  2 Pi  +  Pi-i  =  l~^fdiv  (ui)  ■ 

Considering  the  case  which  has  canopy  interface  between  Ui  and  Ui+ 1,  and  using  Neumann 
boundary  will  give  the  formula 

pAx2 

-Pi  +  Pi-i  =  —^7 -aiv  [Ui )  • 

Introducing  penetration  ratio  7  gives 

p/\orP‘ 

IPi+i  -  (1  +  7)  Pi  +  IPi-i  =  ~^~div  (u*) . 


Update  new  velocity  uses 


u^1  =  u*  - 


2pAx 


(Pi+i-Pi-i) 
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when  came  across  boundary  between  Ui  and  Ui+ 1,  the  new  formula  coupled  with  penetration 
ratio  is 

The  penetration  ratio  7  is  a  dimensionless  parameter  which  controls  the  volume  of  fluid 
going  through  canopy  via  influencing  boundary  condition.  In  the  case  of  zero  porosity, 
the  boundary  conditions  used  for  intermediate  velocity  and  pressure  in  transportation  and 
projection  process  respectively  are 

Ub  'Uintf 

Vpb  =  0 

where  uintf  is  the  velocity  of  interface,  i.e.  parachute  canopy.  This  gives  the  boundary 
condition  of  fluid  at  interface  through  projection  formula 


U*  =  U  H - ~Vp  ^  Ub  =  Uintf 

P 

which  means  the  velocity  of  fluid  at  interface  is  the  velocity  of  the  interface. 

After  coupling  penetration  ratio,  the  boundary  condition  of  intermediate  velocity  is  de¬ 
rived  from  interpolation 


u 


* 

b 


(1-7)  Uintf  +  7 fo+i,  for  Ui 
(1  -  7)  Uintf  +  7 Ui,  for  ui+ 1 


and  the  boundary  condition  for  pressure  is 


Pi+ 1  -  Pi 
Ax 


7V pb ■ 


Substitute  above  equations  into  projection  formula  gives 

Uintf  +  7  iui+ 1  -  Uintf )  +  7 77  Vpb,  for  Ui 
Ub  —  ^  " 

Uintf  +  7  (Ui  ~  Uintf )  +  7^y  Vpb,  for  Ui+l 

The  difference  between  the  fluid  velocity  at  boundary  and  the  boundary  interface  velocity 
is  regarded  as  velocity  due  to  porosity,  which  can  be  represented  as 

i_ \x ^  /\/ 

^ poro  T  ^  ^^b  T  ^ Pb' 

2  p 

Darcy’s  law  describes  the  fluid  through  porous  medium,  which  has  the  form 


K . 


q  = - Vp 

P 
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where  Vp  is  the  pressure  gradient,  k  is  the  intrinsic  permeability  of  medium,  p  is  the  fluid 
viscosity  and  q  is  the  flux.  The  fluid  velocity  is  related  to  flux  q  by  void  fraction  cj) 

-  (i 

V'poro  ,  • 

0 

Equating  the  fluid  velocity  derived  from  numerical  scheme  with  Darcy’s  law  gives 

7  f^p+^fAu 


assume  Ax2  <C  At  gives  simplified  form 


7 


np 

At(j)p 


Figure  10:  Velocity  in  x  direction,  from  left  to  right:  7  =  0,0.25,0.5,0.75 


Figure  11:  Velocity  in  y  direction,  from  left  to  right:  7  =  0,0.25,0.5,0.75 

1.4  Parachute  Simulation  and  Code  Validation 

Our  primary  objective  is  to  conduct  the  fabric-fluid  coupled  simulations  and  validate  it 
against  the  available  parachute  experiments.  Parachutes  differ  in  geometry  and  dimensions 
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Figure  12:  Vorticity,  from  left  to  right:  7  =  0,0.25,0.5,0.75 


of  the  canopies  and  risers.  Among  these  parachutes,  the  T-10  parachute  is  used  by  the  Army 
as  a  personnel  carrier.  This  type  of  parachute  has  a  parabolic  shape  for  its  canopy  with  a 
vent  at  the  top.  The  G-ll  parachute  is  a  cargo  parachute  with  no  vent.  Its  dimensions  could 
vary  and  is  usually  used  as  a  multiple  parachute  system  to  deliver  supporting  equipment. 
We  have  studied  the  single  1/3-scale  G-ll  parachute  inflation.  The  cross  parachute  has  four 
open  side  vents  and  can  be  used  as  a  sports  parachute  and  small  cargo  carrier.  Figure  13 
shows  the  inflation  sequence  of  a  cross  parachute.  These  three  parachutes  produce  different 
patterns  of  airflow  around  the  canopy  and  exert  different  drags  to  the  parachute  system. 

For  comparison  with  the  indoor  vertical  parachute  test,  we  initialized  the  parachute  with 
flat  circular  canopy  of  2.134  m  (7  ft)  nominal  diameter.  We  carried  out  pre-step  running  to 
deform  the  parachute  canopy  and  then  used  the  resulting  geometry  as  the  initial  state  for  the 
drop  test  simulation.  The  initial  skirt  diameter  makes  different  time  graph  for  full  inflation. 
In  this  test,  there  is  a  breathing  motion  due  to  the  over-inflation  of  the  canopy.  It  is  shown 
that  the  canopy  progresses  into  an  over-inflated  shape,  almost  flat,  and  then  contracts  and 
over-inflates  again.  The  breathing  motion  is  an  oscillatory  motion.  The  canopy  appeared  to 
expel  the  excess  air  by  means  of  the  breathing.  In  the  experiment,  the  breathing  motion  was 
also  caused  by  the  constraint  on  the  parachute,  imposed  by  the  guide  wire.  The  breathing 
motion  in  the  simulation  is  smaller  because  there  is  no  vertical  motion  restriction  such  as  a 
guide  wire.  The  terminal  velocity  of  canopy  has  a  good  agreement  with  the  vertical  parachute 
test  results.  The  experimental  data  shows  that  the  descent  speed  rises  rapidly  to  a  peak.  It 
slows  down  while  the  parachute  inflates  and  then  slowly  approaches  a  steady  descent  speed 
of  4.27  m/s  (14  ft/s).  Even  though  the  numerical  solution  cannot  be  compared  with  the 
experimental  data  during  the  initial  period  of  time  in  the  simulation,  we  have  observed  that  it 
reaches  the  terminal  velocity  of  about  3.9  m/s.  The  difference  between  the  terminal  velocity 
in  numerical  simulation  and  the  experimental  data  is  thought  due  to  the  lack  of  porosity  of 
the  canopy  in  the  current  numerical  model.  We  also  simulated  the  1/3-scale  G-ll  parachute 
with  different  canopy  folding  level,  they  all  converge  to  the  terminal  speed  of  3  m/s  in  several 
seconds. 

Parachute  breathing  is  an  oscillatory  motion  caused  by  the  interaction  between  the  fluid 
force  and  the  canopy  at  the  skirt  of  the  parachute.  The  large  difference  between  upper 
and  lower  side  pressure  causes  vibration  in  the  projected  drag  area  of  the  canopy  and  thus 
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the  oscillation  in  descent  speed.  Numerical  study  of  breathing  frequency  is  important  for 
understanding  the  stability  of  the  parachute.  In  our  simulation,  we  initialized  2.134  m  (7  ft) 
flat-circular  parachute  to  compare  with  experimental  data  by  Tutt  et  al.  [40] .  The  breathing 
frequency  in  our  simulation  is  approximately  1.6  -  2.0  Hz  (0.5  s  to  0.6  s  period).  This  is 
in  the  acceptable  range  with  the  experimental  frequency  which  is  2.0  Hz  (0.5  s  period).  If 
we  use  a  parachute  with  larger  size  or  lower  Reynolds  number,  the  breathing  period  will 
increase.  For  example,  the  average  period  of  breathing  for  the  10.7  m  T-10  parachute  is  2.3 
s  in  experiment  [6],  while  it  is  about  2  s  in  our  simulation. 


Figure  13:  Simulation  of  cross  parachute  unfolding  and  inflation.  The  starting  state  of  the 
parachute  is  partially  folded  state.  The  parachute  has  a  small  horizontal  drift  because  the 
folded  state  is  not  perfect  symmetric. 


1.5  Other  Computational  Development 

In  the  last  year,  we  have  accomplished  several  new  tasks  in  the  parachute  simulation  code 
which  enabled  us  for  more  realistic  and  robust  computation  on  the  parachute  problem. 

GPU  Accelerated  Spring  Module  Graphics  Processing  Unit  (GPU)  computing  [13]  is 
to  use  the  GPUs  together  with  CPUs  to  accelerate  a  general-purpose  scientific  and  engineer¬ 
ing  application.  GPU  computing  can  offer  dramatically  enhanced  application  performance 
by  offloading  computation-intensive  portions  of  the  programming  code  to  the  GPU  units, 
while  the  remainder  of  the  code  still  runs  on  the  CPU.  Joint  CPU/GPU  applications  con¬ 
stitute  a  powerful  combination  because  CPUs  consist  of  a  few  cores  optimized  for  serial 
processing,  while  GPUs  consist  of  thousands  of  smaller,  more  efficient  cores  designed  for 
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Figure  14:  The  reinforcement  cables,  or  the  gores,  are  modeled  using  the  same  spring  system, 
but  with  different  spring  constant  along  the  gore  curves.  The  left  plot  shows  fully  opened 
canopy  with  gores.  After  the  inflation,  the  gore  structure  is  clearly  revealed.  The  right  plot 
shows  the  surface  mesh  and  the  details  of  the  gore  curves. 


Figure  15:  Angled  deployment  of  C-9  parachute  with  the  flow.  The  deployment  angle  be¬ 
tween  the  initial  parachute  and  the  direction  of  flow  are  (from  left  to  right)  15°,  30°,  45°, 
60°,  75°,  respectively. 
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massive  parallel  calculations.  Serial  portions  of  the  code  with  logical  comparison  run  on 
the  CPU  while  floating  point  operation  intensive  parallel  portions  of  the  code  run  on  the 
GPU.  The  total  operation  time  recorded  by  the  CPU  clock  and  time  for  GPU  intensive 
computational  part  are  collected  in  Table  3  in  seconds.  From  Table  3  we  can  conclude  that 
the  application  of  the  GPU  has  a  clear  advantage  in  the  computation  of  spring  model.  To 
fully  take  advantage  of  the  GPU  system,  further  optimization  are  needed  to  the  fluid  solver. 
Table  4  shows  the  hardware  structure  of  our  computer. 


Mesh 

Total  Time 

(sec) 

Per  Step 

(sec) 

Without-GPU 

With-GPU 

CPU 

GPU 

80  x  80  x  20 

1345 

152 

1.71 

0.12 

160  x  160  x  40 

8115 

804 

5.02 

0.25 

320  x  320  x  80 

40291 

5349 

11.63 

0.55 

Table  3:  GPU  and  CPU  computing  time  of  three  dimensional  spring  model.  Spring  models 
with  three  different  mesh  sizes  (80  x  80  x  20,  160  x  160  x  40,  320  x  320  x  80)  were  tested 
with  both  pure  CPU  code  and  hybrid  (CPU  and  GPU)  code.  The  hybrid  code  is  15-20  times 
faster  than  the  pure  CPU  code  for  the  computation  intensive  part.  In  general,  the  overall 
speed  of  the  hybrid  code  is  7-10  times  faster  than  the  pure  CPU  code  for  the  spring-mass 
model  part. 


Hardware 

CPU 

Dual  Eight  Core  XEON  E5-2687W,  3.1GHz 

64GB  DDR3 

32KB  x  16  LI  Cache,  256KB  x  16  L2  Cache,  20MB  x 

2  L3  Cache 

GPU 

Dual  Quadro  6000  with  14  multiprocessor,  448  cores, 
1.15Hz 

6GB  global  memory,  64  KB  constant  memory 

48KB  shared  memory  and  32768  registers  per  multi¬ 
processor 

OS 

Fedora  18  with  kernel  3.9.2-200.fcl8.x86_64 

Software 

Compiler 

gcc  version  4.7.2 

CUDA 

CUDA  Toolkit  5.0 

Table  4:  A  Dell  Precision  T7600  Workstation  with  dual  NVIDIA  Quadro  graphics  cards  was 
used  to  set  up  the  test  environment. 

Parachute  gores  are  stitched  by  the  reinforcement  cables.  The  reinforcement  cables  are 
important  structures  on  the  parachute  canopy  surface.  To  accurately  model  the  aerodynamic 
motion  of  a  parachute,  we  need  to  create  a  mathematical  model  which  reveals  the  geometry 
as  well  as  the  material  strength  of  the  canopy  surface  and  the  gores.  The  gore  structures 
in  the  parachute  system  have  important  effect  on  the  stability  of  the  parachute  motion. 
In  our  model,  we  treat  the  reinforcement  cables  as  curves  embedded  in  the  canopy  surface. 
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Young’s  modulus  for  the  fabric  surface  and  the  reinforcement  cable  are  set  at  different  values. 
This  is  realized  by  assigning  different  spring  constants  to  the  surface  mesh  and  to  the  gore 
curves.  The  insertion  of  gores  as  stiffened  interior  curves  in  the  canopy  surface  mesh  is 
demonstrated  by  Figure  1.  Figure  14  shows  that  when  fully  inflated,  the  spring  system 
reveals  both  the  patches  and  the  indentation  of  the  reinforcement  cables.  The  effect  of  the 
gores  on  suppressing  the  vortex  flow  at  the  top  of  the  canopy  is  being  studied.  This  will  be 
published  in  a  new  paper  we  are  currently  working  on. 

Collision  handling  on  parachute  contact  is  through  the  modified  FronTier  library.  The 
computational  module  for  parachute  simulation  is  built  on  the  FronTier  library.  We  have 
used  many  existing  data  structures  and  functionality  for  geometrical  handling.  However,  the 
study  of  parachute,  airfoil  and  other  fluid  structure  based  simulations  requires  major  revisions 
to  the  software  library.  Among  those  needed  are  functions  to  handle  the  non-manifold  surface 
and  three  dimensional  curves  (not  as  boundary  of  a  surface,  such  as  the  string  chords).  The 
original  FronTier  code  was  designed  for  the  study  of  fluid  interface  instability  problems.  For 
these  problems,  collision  surfaces  are  resolved  by  merging  or  bifurcation.  However,  fabric 
surface  can  neither  merge  nor  bifurcate,  therefore  we  need  to  have  functions  which  can 
carefully  deal  with  the  repulsive  contact  and  collision. 

Global  indexing  is  a  new  feature  added  to  the  parachute  code.  The  original  FronTier  had 
to  deal  with  frequent  surface  mesh  optimization  and  topological  reconstructions.  This 
makes  the  parallelization  based  on  global  index  very  difficult.  As  a  result,  the  original 
FronTier  library  relied  on  floating  point  matching  for  parallel  communication.  The  floating 
point  matching  is  not  100  percent  reliable  and  some  more  complicated  algorithms  have  to 
be  implemented  as  reinforcement.  However  for  fabric  surface,  especially  when  spring  model 
is  used,  the  inter-connectivity  and  proximity  of  the  interface  marker  points  are  not  changed. 
Therefore,  global  indexing  is  ideal  for  the  parallel  communication  of  the  surface.  We  have 
added  a  new  parallel  code  for  this  functionality. 

Modularization  is  emphasized  in  our  code  development.  The  parachute  module  is  an 
independent  application  program.  This  new  module  consists  of  four  components:  (1)  the 
initialization  module,  (2)  the  ODE  module  for  the  spring  system,  (3)  the  PDE  module  for 
fluid  dynamics,  and  (4)  the  FronTier  library  for  the  interface  geometry  handling.  A  total  of 
about  15,000  lines  of  code  have  been  written  for  the  first  three  component.  An  additional 
4,500  lines  have  been  added  to  FronTier  to  adapt  the  library  for  the  computation  of  fabric 
surface. 

1.6  Impact  on  Education  and  Training 

The  grant  W911NF1310291  has  partially  supported  two  graduate  students,  Bernard  Moore 
and  Yiyang  Yang.  One  of  them  (Bernard  Moore)  is  an  American  student  from  a  lower  income 
family.  He  is  a  highly  motivated  student  and  is  interested  in  the  US  defense  research  after 
graduation.  He  participated  in  the  US  Air  Force  SFFP  program  and  is  highly  motivated 
by  our  project  as  well  as  the  service  to  the  military.  Another  former  student  from  our 
research  group,  Eric  Rathnayke  has  joined  Military  Accessions  Vital  to  the  National  Interest 
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(MAVNI)  program  this  year.  The  numerical  techniques  we  studied  so  far  have  had  an  impact 
on  education  and  training. 

Yiyang  Yang  will  soon  graduate  with  the  Ph.  D.  degree  in  computational  and  applied 
mathematics.  She  has  already  been  offered  as  a  software  engineer  by  Goldman  Sachs.  An¬ 
other  student,  Qiangqiang  Shi,  not  supported  by  this  grant  but  played  an  important  role 
in  this  project,  has  already  received  his  Ph.  D.  degree  and  joined  Bloomberg  as  a  software 
engineer. 

1.7  Communication  and  Outreach  of  the  Project 

Our  project  has  been  constantly  advised  by  Army  scientist  Dr.  Richard  Charles  with  his 
expertise  in  modelling,  experiment  and  field  experience  on  the  parachute  delivery  system. 
One  joint  paper  has  been  submitted  to  the  Journal  of  Fluid  and  Structure  which  is  still 
under  review.  We  anticipate  a  second  paper  with  him  in  the  coming  months.  The  PI  has 
also  reached  Air  Force  scientist  Dr.  Daniel  Reasor  in  the  Air  Force  Summer  Fellowship 
program.  Other  collaborations  include  a  joint  research  project  with  Drs.  Robert  McGraw 
and  Yangang  Liu  at  Brookhaven  National  Lab  on  climate  modelling  and  with  Dr.  Valmor 
de  Almeida  at  Oak  Ridge  National  Lab  on  the  study  of  phase  transition  problems  in  nuclear 
reactor. 
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